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1.0 Introduction 



In a severe thermal environment, materials generally ex- 
perience a degraduation of stiffness and strength, and as a 
result there is a diminishment in the load carrying capacity 
of the material. In order to design a structure which will 
perform safely in a thermal environment, the structural engi- 
neer must have an analysis program which will predict the be- 
havior of the structure. This report documents an analysis 
for the behavior of a fiber reinforced composite plate sub- 
jected to loads in a thermal environment. 

2 . 0 Description of Fiber Reinforced Composite Plates 

As a result of imperfections, such as dislocations and 
surface scratches, the strength of bulk material is orders of 
magnitude less than the strength of the material in crystal 
or whisker form. The cost of producing crystals is so great 
as to prohibit production except on a limited basis. On the 
other hand, fibers of a material can be produced at a cost 
which permits their application to structural systems. 

Depending upon it's diameter, a fiber of a material has 
a strength intermediate to the bulk and crystal strengths of 
the material. For example, the strength of steel fiber is 
approximately 4.14 MPa (600,000 psi), compared to a strength 
of 1.04 MPa (150,000 psi) for steel in it's common bulk form. 
Thus the use of steel fiber in place of common steel results 
in a four-fold reduction in weight. Actually the use of fiber 
requires their embedment in a matrix of another material in 
order to achieve an useful structural form. However even with 
embedment, fiber reinforced materials offer a significant in- 
crease in the strength to weight ratio over common materials. 
This is the primary reason for using composite materials. 

One of the most common useful structural component is 
the plate element. A fiber reinforced composite plate (or 
simply, a laminate) consists of layers of lamina (or plies) 
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stacked one upon another. A lamina itself is comprised of uni- 
directional fibers embedded in a binding material, called the 
matrix. Generally, the binding matrix material is a much weak- 
er and more ductile material than the fiber reinforcement em- 
bedded within it. Some common fiber/matrix combinations are 
graphite/epoxy , glass/epoxy, and boron/aluminum. Some idea of 
the dimensions involved can be obtained from the dimensions 
of a typical graphite/epoxy lamina; 

i) fiber diameter is approximately .008 mm. (.0003 in.] 

ii). lamina thickness is approximately .125 mm. (.005 in.] 

iii] fiber volume to total volume (called the fiber 
fraction) is about 0.5 

For these dimensions, a 12.5 mm. (0.5 in.) thick laminate 

will consist of 100 plies (lamina) stacked one upon another. 

As a result of unidirectional fibers embedded within it, 

a lamina is not isotropic. The stiffness of a lamina along it's 

fiber direction (called the longitudinal stiffness) is much 

greater than the stiffness of the lamina in the direction 

orthogonal to the fiber (called the transverse stiffness). 

For example, in the case of a typical graphite/ epoxy lamina, 

the longitudinal stiffness (E or E^) and the transverse 

stiffness (E or E r ) are 181 GPa (26x10^ psi) and 10.3 GPa 
6 > 1 

(1.5x10 psi), respectively. 

Stacking lamina of different orientation and thick- 
ness upon one another to form a composite plate results in an 
orthotropic body. One of the most attractive features of a 
composite plate is the capacity for design of material prop- 
erties, an option not previously available to structural 
engineers . 

3. 0 Lamina and Laminate Coordinate Systems 

In dealing with fiber reinforced composite plates, it is 
necessary to employ two types of coordinate systems. The first 
type of coordinate system is associated with a lamina. This 
coordinate system is referred to as the natural coordinate 
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system, or the local coordinate system, or the lamina (or ply) 
coordinate system. Here we assign the cartesian pair of x,y 
axes to the plane of the lamina, with the x-axis along the 
fiber direction of the lamina. Since adjacent lamina have 
different fiber orientation with respect to one another, there 
are as many natural coordinate systems as there are plies. The 
i-th lamina (ply) has the coordinate system (x,y)^ or (x^,y^). 
This is shown in Fig.l below. 




Figure 1. Lamina Layup to Form a Laminate 

The second type of coordinate system associated with a 
composite plate is the coordinate system which ties all of 
the lamina coordinate systems to a single pair of axes, 
called the laminate or global coordinate system, and denoted 
by (x,y). These axes lie in a plane parallel to the lamina 
planes, and are generally located at the midplane of the la- 
minate, as shown in Fig. 2 for a 3 ply laminate. 

Each of the local (lamina) coordinate systems are rela- • 
ted to the global (laminate) coordinate system by the angle 
4>^, where <fu is the angle between the x^ axis of the i-th 
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lamina and the x axis of the laminate. In this way, the 
orientation layup of a n-ply laminate can be denoted by the 
notation [<j> 1 , 4*2 , 4*3 4> n ] . There is no requirement that 




Figure 2. Lamina and Laminate Coordinate Systems 

the lamina be of equal thickness, although this is a common 
construction. 

4 . 0 Stress-Strain Relations in Natural Coordinates for a 
Lamina in Plane Stress 

In the case of plates subjected to inplane and trans- 
verse loads, the result is the plane stress state. This means 
that for a lamina in the x,y plane, interlaminar stresses <? 2Z , 

a zy’ an< ^ a zx are negligible compared to inplane stresses a , 

AA 

a , and o , and are therefore ignored. It should be noted 
in passing that interlaminar stresses can be very large in 
local regions of boundaries. (See Ref.l for details). In this 
report interlaminar stresses are not considered. 
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In terms of it's own natural coordinate system, the 
stress-strain relations for the i-th lamina, in tensor nota- 
tion, are 



j i krs e 



rs" e jk 0) i j.k,r,s=x 1 ,y i (1) 



where a.,, and e are the second order stress and strain 
j k rs 

tensors, 0^ is the second order tensor of thermal coeffi- 
cients, C., is the fourth order tensor of stiffness 
J s 

coefficients, and 0 is the difference in temperature from 
the stress free temperature. That is, if T is the stress free 
(curing) temperature, and is the temperature of the i-th 
lamina, then 

9 i = (T i - T) (2) 

Throughout this report, the Einstein summation convention 
is in effect for all tensor equations. 

Alternatively, Eq.(l) may be written for strains in 
terms of stress, as 

Ce jk ) i' (S jkrs < ’rs* a jk 0) l j.k.r.s-x^ (3) 

where Sjj crs is the fourth order tensor of compliance coeff- 
icients, and a j ^ is the tensor of thermal expansion 
coefficients. We note that properties of a lamina, such as 
stiffnesses, compliances, thermal coefficients, and thermal 
expansion coefficients, are functions of temperature. How- 
ever over narrow ranges of temperature, properties can be 
treated as constant without introducing significant error 
into the analysis. 

As a matter of convenience, hereafter we proceed 
using contracted "vector" notation. It should always be 
kept in mind that stresses, strains, thermal coefficients, 
thermal expansion coefficients, stiffnesses, and compliances 
are, in fact, tensors, not vectors, and therefore obey the 
appropriate tensor transformation laws when transformations 
occur. The second order tensors ,e ij ' an< ^ a ij^ * n 

the contracted notation have a single subscript in accord- 
ance with Table 1, and the fourth order tensors and 

Cfjki) carry the double subscript notation shown in Table 2. 
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Tensor 


Contracted 


Subscripts 


Notation 


XX 
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yy 


2 


xy and yx 
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Tensor 

Subscripts 


Contracted 

Notation 


xxxx 




11 


xxyy 




12 


xxxy 


and xxyx 


16 


yyxx 




21 


yyyy 




22 


yyxy 


and yyyx 


26 


xyxx 


and yxxx 


61 


xyyy 


and yxyy 


62 


xyxy 


and yxyx 


66 



Table 1. Contracted Notation 

for 2 n< ^ order Tensors 



Table 2. Contracted Notation for 
"th. 

4 order Tensors 



For example, using the symmetry conditions, e *e , 

xy yx 

C xxxy = Scxyx ancl s0 on > anci defining the engineering shear 
strain Y *(e + el , equations (1) become 

Ay Ay y a 



[a l = C ll e l + C 12 e 2 + C 16 y 6 " B l 9 ^i 

[ a 2 = C 21 e l + C 22 e 2 + C 26 Y 6 ' g 2 0 ] i 

t a 6 = C 61 e l + C 62 e 2 + C 66 y 6 ' S 6 9 3 f 



or, in vector notation 

o.= ( Ce - 39 ) (5) 

X ~ ~ ~ X 

where the i subscript denotes the equation is for the i-th 
lamina. In equation (5), and hereafter, a single tilda under 
a quantity denotes the quantity is treated as a 3x1 vector, 
a double tilda under a quantity means the quantity is treated 
as a 3x3 matrix, and a quantity without a tilda beneath it 
denotes that the quantity is a scalar. 

As a result of the natural coordinates being along and 
orthogonal to the fibers of a lamina, the natural coordinates 
are the principal axes of the lamina. Thus the normal and 
shear effects uncouple, and Eq. (5) takes on the special form. 



6 



0 





66 



2 

Y 6 



\ 







i i i 



( 6 ) 



That is, in the natural coordinates of a lamina, 

=0, and $6=0. Furthermore, as the xy plane is a plane of 
symmetry, C]_ 2 =E 2 i> an< ^ therefore C is a symmetric matrix. 
Thus a lamina has four independent stiffness coefficients 
(Cn , Cf 2 , C 22 > and ^66^’ an< ^ two independent thermal coeffi- 
cients C 3 2 . , an< ^ ^ 2 ^ ‘ Such materials are called orthotropic 
materials. It is important to note that the special form of 
C in Eq.(6), associated with uncoupled normal and shear be- 
havior, is valid only in natural coordinates of a lamina. 

The notation used most frequently in the composites 
literature is § for compliance, and 0 rather than C for 

ft 

stiffness. Hereafter we shall use S and Q. In contracted 
notation, Eq's.(l) and (3) are, 



o. = (Qe -ee^ (7) 

and 

Si = (So +a0) i (8) 

Premultiplication of Eq.(7) by Q ^ and solving for strain 
e gives 

Si = CQ _1 a + Q' 1 ee) i ( 9 ) 

A comparison of equations (8) and (9) shows the relations 
between Q and S, and § and a are 

Q = S' 1 (10) 

and a = Q ‘*'8 . (11) 

The stiffness coefficients Q. v are related to the common 

1 * 

engineering coefficients E x > E y> v X y> and ^ X y» E y Eq's.(12). 

The engineering coefficients are obtained by performing the 
simple tension and shear tests on a lamina as depicted in 
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Figure 3. Pulling a lamina along the fiber direction (Fig. 3a) 
gives the longitudinal stiffness, E . The tension test in the 

A 

orthogonal direction (Fig. 3b) gives the transverse stiffness 

E . The longitudinal Poisson's ratio v is the ratio of e 
y 6 xy x 

to e^ obtained in the longitudinal tension test. The long- 
itudinal shear modulus G is obtained from the shear test 

xy 

shown in Fig. 3c. The relations between S., and engineering 

j K 

coefficients are, 



S n -(1/E x ) 

s i2 =s 2r-(V E x )= -tv /Ey) tl2) 

S 22 =U/E y ) 

S 66=t ly V 

Note that the second of these relations defines a transverse 

Poisson's ratio, v . The x and y subscripts denoting longi- 

yx 

tudinal and transverse properties are quite often replaced by 
T and L subscripts, i.e., E^ denotes the longitudinal stiff- 
ness of a lamina, E^. is the transverse stiffness, and so on. 
Using Eq. (10) , the relations between stiffness coefficients 
Q— and the engineering constants are found to be. 



Qu-cy^-'W” 



^12~^21 ^ v xy^x^ v xy v yx^ ^ v yx^y^^ v xy v yx^ 



Q„=(E /(1-v v 1 



: 2 2 v y 



xy yx' 




(13) 



Figure 3. Experimental Tests for Engineering 
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5 . 0 Stress-Strain Relations for a Lamina in Plane Stress 

with Respect to Laminate Axes 

In order to develope the equations for a laminate, it is 
necessary to relate each lamina comprising the laminate to a 
common coordinate system, the laminate axes. This is accomp- 
lished by coordinate transformation from lamina axes to lami- 
nate axes. Each lamina coordinate system (x,y,z)^ is oriented 
to the laminate axes (x,y,z), where z. and z coincide, by an 
angle <}k. Recall from Fig. 2, that is the angle between the 
x^ axis of the i-th lamina and the x axis of the laminate. 

Hereafter all quantities with a bar over the quantity 
symbol means that the quantity is with respect to laminate 
axes, and quantities without a bar over them means the quant- 
ities are with respect to lamina axes. Thus the stresses, 
strains, thermal coefficients, and stiffnesses with respect 
to the laminate coordinate system are a,e,g, and Q. Using the 
tensor transformation laws for the second order tensors. 



where 1 is the cosine of the angle between the r-th lamina 
axis and the s-th laminate axis, gives the expressions for 
lamina quantities with respect to laminate axes. Denoting 
the transformation between coordinate systems by the matrix 
T, Eq's.(14) may be written, 




(14) 






(15) 



Si » (Tg)i 
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Substitution of Eq's.(15) into Eq's.(7) gives 



T . a . = Q.T.e. - T • f3 • 0. 

Sl-l *»1~1 1 

Premultiplication through by the inverse of gives, 

a. = TT 1 Q.T. e. - B.0. 

■” x »i 1~1 K 1 1 

Denoting the matrix by Q^, where is the 3x3 stiff- 

ness matrix of the i-th lamina with respect to laminate axes, 
we obtain the stress-strain relations for a lamina in terms 
of laminate axes as, 

o .= Q . e - 3.0. 

*1 5lSi Cl l 

Recalling Eq.(ll), yields the more useful expression (since 
the coefficients of thermal expansion are generally known) , 



2i* Si 9 ! 5 



(16) 



where the i subscript refers to the i-th lamina. Explicitely 
Eq' s . (16) are. 



(17) 



V 






$12 


*16 




® i " 




V 


*2 




§21 


^22 


*26 




e 2 


- 0. 
1 


“2 


, a 6. 


i 


.^61 


^62 


^ 66 _ 


i 


. Y 6. 


i 





where, upon using c^= cos<f>^ and s^= sin<f>^, we obtain, 

CQllV [Q n c 4 + 2(Q 12 + 2Q 66 )s 2 c 2 + Q 22 s 4 ] i 

^12^ i = ^ ( - ( ^ll + ^22" 4( ^66^ S C + < ^12 ( - s +C ^i 
^22^i = ^ll s +2 ^12 + 2( ^66 ^ s c + ^22 c ^ i 
^16^i = f ( ‘ ( ^ll'^12" 2( ^66-' SC + ^12'^22 + 2( ^66- ) s i 
^26-*i = f^ll'^12' 2 ^66^ s c + ^12"^22 + 2( ^66-* SC ^ i 

CQ 6 eV t(Qn + Q 2 2- 2 Qi2- 2 Q66) s2c2 + Q 66 ^ 4+c4 ^i 



(18) 
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We note from Eq’s.(18) that the stiffness matrix of a 
lamina with respect to the laminate coordinate system is a 
full 3x3 matrix. However, there are still only four inde- 
pendent coefficients; any two of the coefficients may 

be written as functions of the remaining four coefficients. 
The Q matrix is symmetric, that is, Q-^ = $21’ Qjg = $ 61 > 

and ^26 = $62* 

The fact that Q is, in general, full, means that coupling 
between normal and shear behavior occurs. In contrast, for 
isotropic materials, normal and shear behaviors uncouple. 

For cross ply laminates, that is laminates with orthogonal 
plies, subjected to inplane forces and bending moments, 
normal and shear behaviors uncouple. For angle ply laminates, 
that is laminates with some non-orthogonal plies, coupling 
between normal and shear behaviors occurs. This difference 
in behavior between cross ply laminates and angle ply lami- 
nates is shown by an illustrative example in section 11. 



6 . 0 Strain-Displacement Relations for a Laminate 

In accordance with common thin plate theory, the as- 
sumption that planes before loading remain planes after load- 
ing (Kirchhoff-Love theory) is taken as the basis of deform- 
ation of a laminate. In this case, the displacement field is 
given by, 



u = u - z(3w/3x 



V = 



v Q - z(8w/3y) 



0 ) 



(19) 



where u and v are the displacements in the x and y directions, 
u q and v Q are the midplane displacements, w is the displace- 
ment in the z direction, and z is the distance from the mid- 
plane to the point of interest. 

Using the linear strain-displacement relations, 



e = ( 9u/3x) e = (3v/3y) 

a y 

Y X y = ( 3u/ 3y) +,(3v/3x) 



( 20 ) 
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equations (19) become. 



e l = ( 3 u c/ 3x )" zf? 2 ™/ 9 *, 2 ) 

e 2 = (3v o /3y)- z(3?w/3y 2 ) l (21) 

Yg = (3u Q /3y + 3 v q /3x) - 2 z (c^w/ 3x3y) 

Equations (21) may be written in vector form as, 

e = e Q + zk (22) 

where e Q is the strain vector associated with midplane exten- 
sional displacements given by 



(Vo 


3u q /3x \ 




( Vo ‘ 


3v Q /3y k 


(23) 


(Vo 


3u Q /3y + 3v q /3x ) 





and ic is the curvature vector associated with bending of the 
midplane, given by 



(k x ) 


- S^/^x N 




(< 2 ) 


- 3 2 w/3 2 y » 


(24) 


(k 6 ) 


-23 2 w/3x3y / 





Since the midplane strains and curvatures belong to the lami- 
nate, they do not change from one lamina to another, and 
therefore the bar over these quantities is not required and 
henceforth will be dispensed with. 



7 . 0 Force-Displacement Relations for a Laminate 

Equation (22) shows that the strain at a point is the 
sum of two components, (i) midplane extensional strain, and 
(ii) midplane curvature. Figure 4. depicts the displacement, 
strain, modulus, and stress field with respect to the x,z 
plane at an arbitrary y location. The displacement and 
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strain fields are shown together in Fig. 4b, while in fact, 
one is a scalar times the other. In accordance with the 
Kirchhoff-Love assumption, the displacement field is linear 
across the laminate thickness. As a result of the jump dis- 
continuity of modulus from lamina to lamina (Fig. 4c), the 
stress field also has a jump discontinuity from lamina to 
lamina (Fig . 4d) . 




V Y 



Laminate 

Field 




Displacement (and 
Strain) Field 




(a) (b) (c) 



(d) 



Figure 4. Displacement, Strain, Modulus and Stress Distribu - 
tions 



The laminate forces N and moments M associated with a la- 
minate are obtained by integrations of the stress field over 
the laminate thickness. The extensional laminate forces are 
given by, 




( 25 ) 
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and the laminate moments are given by, 



- *1 
M- 


1 

/•t 


" _ - 

0, 


1 


1 2 


1 


M 2 


-J 


5 2 


**6 

L d 


•/ 

- t 
* 
L 


u & - 



zdz 



y -Jl 



idz (26) 



where t is the plate thickness. The laminate forces and mo- 
ments are shown in Figure 5. Components and ^ are the 
normal components of force, and is the shear component. M^ 

and M 2 are bending moments, and M^ is the twisting moment. 




Figure 5. Laminate Forces and Moments 



Since the stress field across the laminate thickness is 

the sum of the lamina stresses, 

n 

5-to. 

i=f 1 

equations (25) and (26) become 

.z. 



and 



* -fc/i 

5.-1 

- ^ C- z± 

M =ZJ o.z. 



dz 



(27) 



(28) 



z i- 1 
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where z^_^ and z^ are the upper and lower limits of integra- 
tion for the i-th lamina, and n is the number of lamina com- 
prising the laminate. 

Substituting the expressions for lamina stresses in 
terms of strains, Eq's.(16), into the previous Eq's.(27) and 
(28) gives the laminate forces and moments in terms of strains 
as , 





N = 


.is l 9itSi- 


a i 9 i )dz 


(29) 


and 


M ■ 


1=1 H-i 

.2; Zi 3i(§i- 

i-izi-i 


a^9^) zdz 


(30) 



The strain in the i-th lamina is given by Eq.(22), 



e. = e„ + z k (31) 

~i ~o 

where e Q and k are the midplane extensional strain and mid- 
plane curvature vectors respectively, and are constant from 
one lamina to another. Substituting Eq.(31) into Eq's.(29) 
and (30) gives, 

z . 

N = § / 1 Q ± ( e Q + zk - a i 9 i )dz 

i=lz i-*l 

nZ-_ - -- - _ _ 

M = z iQ-Ce + zk - a.9.)zdz 
•_,j ~1 ~0 - -1 l' 

1_lz i-"l 

Since e Q and k do not vary with z, and and 9^ are con- 

stant for a lamina, the previous equations become, 

n z j, _ z i_ 

N = z J“ z + K J zdz l 

1=1 ~ ' z i-i ' z i-i 

n _ Zj__ _ _ _ 

M = I Qi [ ( e o -c *i 0 i) / zdz + k J zdz] 

1=1=1 ' Z i-1 ~ z i-i 

Upon integration, we obtain, 

N « Z QiC® o "“i 0 i^*i-*i-l^ + 1 ( 32 ) 

i=l ~ i=l ~ 2 
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and 



M = 



n 



J'i * i= 2 yj( yvi) 



Defining the matrices, 



and vectors , 



a = ? sPvy-P 

i=l - 


(34) 


B = r Q- (z?-z? .. )/2 
. “ 5i i l-l 
1=1 


(35) 


D = ? Qi(z?-z3 )/6 
~ i=l ~ 


(36) 




(37) 


“ Z Q.a.0.(z?-z? ,)/2 

i=l - 1 - 1 1 1 1-1 


(38) 


(33) become, 




= Ae + Bi< - N t 


(39) 


= Be Q + Dg - Mj 


(40) 



Equations (39) and (40) are the "force- displacement" rela- 
tions for a laminate. 

The and Hj, vectors are the so called thermal force 
vector and thermal moment vector respectively, as they result 
from a thermal environment. If the temperature of each lami- 
na is equal to the stress free temperature (i.e., the curing 
temperature) , then 9^ is zero for each lamina, and the therm- 
al vectors are equal to the zero vector. In this case, Eq's. 
(39) and (40) reduce to the force- displacement relations for 
a laminate not subjected to thermal effects. 

The A matrix relates inplane forces to inplane strains, 
and the D matrix relates moments to curvature. The B matrix 
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couples inplane forces to curvature, as well as moments to 
inplane strains. In the special, albeit common, case of sym- 
metric laminates, the contribution to 6 from a lamina above 
the midplane is cancelled by the contribution from a counter- 
part lamina below the midplane, and as a result, the B matrix 
is the zero matrix, i.e., 6=0. Thus for symmetric laminates, 



o" 


n t 


(41) 


Dk - 


m t 


(42) 



and therefore the extensional and bending behaviors are un- 
coupled. If the laminate is symmetric and the temperature 
distribution is symmetric with respect to the midplane, then 
from Eq.(38), the thermal moment vector is equal to the zero 
vector and Eq. (42) reduces to 

M = Dk ' (43) 



8 . 0 Stress-Force Relations for a Lamina in a Symmetric 
Laminate 

In order to determine whether a ply has failed, it is 
necessary to determine the stresses in each ply. The devel- 
opement of the stress-force relations for a lamina, which 
follows, is restricted to symmetric laminates. This is not a 
severe restriction since symmetric laminates are in most 
common use. Solving Eq's.(41) and (42) for strains and curv- 
atures gives, 



e = A _1 (N + N t ) 
- o ~ ~ ~ I 

and 

k = D"' l (M + m t ) 



(44) 

(45) 



The lamina strains are related to the midplane strains 
and curvatures of the laminate by Eq.(31), 



e . = e + zk 
~ 1 ~o 



( 31 ) 
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Substitution of Eq’s.(44) and (45) into Eq.(31) gives, 

§i = A~ 1 (N + N t ) + z i p' 1 (M + Rj.) (46) 

where z^ is the distance from the laminate midplane to the 
middle of the i-th lamina, that is, 



z . 

l 




z 



i-1 



)/2 



(47) 



Lamina stresses with respect to laminate axes are related to 
lamina strains in laminate axes by Eq.(16) as. 



a. = Q.(e. - a.0.) 

~i ->i ~i ~i i' 



(16) 



Finally then, substitution of Eq.(46) into Eq.(16) gives the 
lamina stresses, with respect to laminate axes, in terms of 
laminate inplane forces N, laminate moments M, thermal forces 
N-p, and thermal moments as, 



Of= Q i [A’ 1 (N + N t ) + ^D'^M + M t ) -o i 0 i ] 



(48) 



Equation (48) takes on special forms for particular cases; a 
few of which are presented below. 

a) Temperature free effects occur when the plate is at a 
uniform temperature equal to the curing temperature of 
laminate. In this case, 0. =0, and therefore N„ and M™ 
are equal to the zero vector, and Eq.(48) reduces to, 

a i = Q^A" 1 ^ + z i D _1 M) (49) 

In the absence of lateral loads on the laminate, M = 0, 
and Eq.(49) becomes, 

o = QA -1 N (50) 

b) Absence of moments occurs when the laminate is not sub- 
jected to lateral loads. In this case, M = 0, and 
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equation (48) reduces to, 




( 51 ) 



If the temperature field is symmetric with respect to 
the midplane of the laminate, then Rp is the zero 
vector and Eq.(51) becomes, 




(52) 



c) Thermal effects only occur when the laminate is not 
subjected to inplane and lateral loads. In this case 
N = M = 0, and Eq.(48) reduces to 



If the temperature field is symmetric with respect to 
the laminate midplane, then Rp = 0, and Eq.(53) gives, 



9 . 0 Failure Criterion 

When a composite plate is subjected to a load and/or 
thermal environment, stresses develope in each lamina in ac- 
cordance with Eq.(48). Depending upon the magnitude of the 
lamina stresses, one or more lamina may fail. When a lamina 
fails, there are two effects. First, if we assume brittle 
failure, then a lamina will relinquish it's load carrying 
capacity upon failure, and there will be a redistribution 
of load to the remaining non failed lamina. Secondly, when a 
lamina fails, the stiffness matrices A and D change in ac- 
cordance with Eq's.(34) and (36), and thus the redistribu- 
tion of stresses in the non failed lamina is not simply a 
change in magnitude, but rather a change in the stress state. 

The hypothesis of a brittle failure of a lamina is an 




(53) 




(54) 
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idealization of the actual behavior. Actually, lamina have 
some ductility and therefore continue to carry some fraction 
of load, and contribute some stiffness to the laminate upon 
failure. Here we assume a failed lamina serves no useful pur- 
pose after it's failure. 

As will be demonstrated by illustrative examples in 
section 11, failure of a ply does not necessarily mean that 
the laminate has failed. There we will show that a cross ply 
laminate subjected to uniaxial loading does not fail upon 
failure of the first ply, and therefore the magnitude of the 
load can be increased beyond first ply failure. We also show 
that under uniaxial loading, angle ply laminates fail when 
the first ply fails. 

There is no way to establish the validity of any failure 
criterion on theoretical grounds. All that can be done is to 
show that a failure criterion gives analytical results which 
agree, for a wide class of problems, with experimental re- 
sults. There are several failure theories which give good re- 
sults and are in common use. The most general of the macro- 
scopic failure criteria, which includes many of the other 
criteria as special cases, is the Tsai-Wu criterion. Ref. (2). 
It is the failure criterion adopted in this work. 

According to the Tsai-Wu criterion, a lamina fails when 
the general quadratic equation (in tensor form) , 



F ij°i a j + F i a i = 1 j = 1 > 2 , , , 6 (55) 

is satisfied. In Eq.(55), and F^. are strength tensors of 
second and fourth order, and the equation is in contracted 
notation. In the case of plane stress, the stress state for 
composite plates, Eq. (55) in terms of natural coordinates of 
a lamina, becomes 



F ll a l + 2F 12°l a 2 



F 22 a 2 + F 66 a 6 + F l a l 



+ F 2 a 2 - 1 



( 56 ) 



Equation (56) is applied to each lamina. If the left hand 
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side of Eq. (56) is less than unity, then the lamina has not 
failed. If the left hand side is greater than unity, then the 
lamina has failed. Failure of a lamina is imminent (impend- 
ing) when the left hand side is equal to unity, i.e., when 
Eq. (56) is satisfied identically. 

In order to use Eq.(56) directly, that is without modi- 
fication, it is necessary to use lamina stresses with respect 
to lamina axes. The use of lamina stresses with respect to 
laminate axes would require the transformation of the second 
and fourth order tensors, and F^j , to laminate axes. Al- 
though this could be done, it is more expeditious to trans- 
form the lamina stresses with respect to laminate axes, 
given by Eq's.(48), into lamina stresses with respect to la- 
mina axes. The transformation is accomplished by the use of 
Eq.(14). Thus, we obtain, 



(cPi . (oPicf ♦ OPjS? * 2(5 6 ) i c i s i 
(a 2 Ji = (Spp? ♦ (o 2 J iC ? - 2(5 6 ) i c i s i 



(57) 



(a 6 )i -CopcppCip.CjS. + CSppcI-sf) 



where c. = cos<f>., s. = sin<J>., and i denotes the lamina. 
i i’ i Y i* 

Five of the six strength parameters of Eq.(56) may be 
determined by a series of simple experiments as follows. A 
tension test of a lamina along it's fiber direction gives 
the longitudinal tension strength, X. A compression test 
along the fiber axis gives the longitudinal compressive 
strength, X* . Tension and compressive tests in the transverse 
direction (i.e., orthogonal to the fiber), gives the trans- 
verse tension and compressive strengths, Y and Y', respect- 
ively. A shear test of the lamina gives it's shear strength, 
S. The strength parameters in Eq.(56) are related to the 
experimentally determined strengths as follows. Applying 
Eq. (56) to each of the five stress states associated with 
the five tests gives five simultaneous equations which may 
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be solved for F n > F 2 2 ’ F 66 ’ F 1 ’ an< ^ ^2 For exam pl e > Eq. (56) 
for the longitudinal tension test becomes, 

F 11 X2 + F i x “ 1 

and for the longitudinal compression test, Eq.(56) is, 

F 11 X,2 ‘ F i x ' = 1 

Solving these equations for and gives the first two 

of Eq's.(58). The remaining strength parameters are obtained 
in a similar way. 



F n = ( 1/XX ' ) 

F x = (1/X) - (1/X') 

F 22 = ( 1/YY 1 ) (58) 

F 2 = (1/Y) - ( 1/ Y * ) 

F 66 = ^/ S2 ^ 

The remaining strength parameter F 12 is obtained on 
mathematical grounds. In order for the failure to be a closed 
curve, the quadratic equation must be an ellipse. This leads 
to the condition (see Ref. (3) for details), 

-1 < FJ 2 /F u F 22 < 1 

If the Tsai-Wu criterion is required to be a generalization 
of the von Mises criterion, then ( F i2/ F 11 F 22^ must ec l ua l 
to -0.5. This condition leads to, 

F 12 = * 5 *- F ll F 22'* ^ 59 -* 

Some typical values of lamina strengths are presented in 
Appendix A. 

10.0 Strength Ratio and Failure Analysis 

As discussed in the previous section, the Tsai-Wu cri- 
terion of Eq.(56) can be used to determine whether or not a 
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ply has failed by evaluation of the left hand side of the 
equation. By introducing a parameter, called the strength 
ratio, the criterion can be used to determine when a ply 
fails, the sequence of ply failure, and the load for laminate 
failure. 

Given the temperature field of a laminate, T^, and the 
relative magnitudes of the components of the laminate "load" 
vectors, N and M, we seek to determine the allowable load 
vectors, N and M , which will bring a lamina to imminent 
failure. Since the laminate stiffness matrices, A and D, 
change with every ply failure, a cyclic (iterative) pro- 
cedure is required to trace the sequence of ply failure. 

The m-th cycle determines the m-th ply which will fail. 
During the m-th cycle, the strength ratio of the i-th ply is 
defined as, 



R i(m) 



= N 




1 



(60) 



This ratio is the same for all components of N and M, that 

is, N. /N. = M. /M. , for j-1,2,6. The allowable stresses in 
’ ja / j ja' j ’ J ’ ’ 

a lamina during the m-th cycle, with respect to laminate axes, 
are given by Eq's.(48), 



^■?ia^ (m) 



2i [ Fw ( 5aW + 5T )+i i5(m)«a(m)* 



- a. 6. ] 



(48) 



where the 'a 1 subscript denotes allowable, the i subscript de- 
notes the i-th lamina, the m subscript denotes the m-th cycle, 
and A r •> and i) r -> are the stiffness matrices of the m-th 
cycle laminate. In Eq. (48) , i takes on the values of only the 
non failed lamina. The m-th cycle laminate is comprised of 
the lamina which have not failed during the previous (m-1) 
cycles . 



It is convenient to partition the lamina stresses into 
two parts; one part due to the N and M loads, and another 
part due to the N^. and WLp thermal vectors. Denoting the 
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stresses due to loads by 0 .^, and the stresses due to temper- 
ature by we have, 



G ia(m) 0 iL(m) + 0 iT(m) 



(61) 



where the load stresses 2iL(m)’ an< ^ thermal stresses 
are given by, 

b.,, -> = Q . [ A^ 1 -, -N r , + z.D*A *M , .] (62) 

~iL(m) ~i ~(m) ~a(m) i~ (m) ~.a(m) J v 



~iT (m) ' Si t J(i)’?T + z i?(m)'?T ;i 



(63) 



Using Eq's.(57), the above lamina stresses with respect to 
laminate axes are transformed into lamina stresses with re- 
spect to lamina axes for use in Eq.(56). Substituting lamina 
stresses with respect to lamina axes into the Tsai-Wu cri- 
terion of Eq.(56) gives a quadratic equation for the i-th 
lamina and m-th cycle, 



R 2 

a i(m) R i(m) 



+ b . , % R. , + c . =0 

i(m) i(m) i(m) 



(64) 



where the coefficients are. 



a i(m)“ [ F ll a lL + 2F 12 a lL a 2L + F 22 a 2L 



+ F 66 G 6L J i(m) 



(65) 



b i(m) = ^ 2( ' F ll a lL 0 lT + F 12 a lL G 2T + F 12 G 2L G 1T 

+ F 22 0 2L a 2T + F 66 a 6L a 6T^ ( - 66 ' ) 

+ F 1 0 1L + F 2 G 2lJ i(m) 

C i(m) = ^ F 11 0 1T +2F 12 0 1T 0 2T + F 22 G 2T 

(67) 

+ F 66 g 6T + F 1 g 1T + F 2 g 2T^ i(m) ' 1 
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The strength ratios for each of the lamina compris 

ing the m-th cycle laminate are obtained as the roots of Eq. 
(64) . The lamina with the smallest absolute strength ratio is 
the lamina which will fail during the m-th cycle. Thus the 
load factor to bring the m-th ply to failure is, 

R W * Mi„( Ab sCR i(nl ))l (68) 

Laminate failure occurs when all lamina have failed. The 
strength ratio (or load factor) for the laminate to fail, de- 
noted by R u , is the maximum of the load factors, that is 

R u = Max[R* (m) ] (69) 



If the maximum strength ratio is obtained for m<n, where n is 
the number of plies comprising the laminate, then the lamin- 
ate fails after m plies have failed. The allowable loads that 
the laminate can sustain without failure, for a given temp- 
erature distribution, are obtained from Eq.(60) as, 




R U N 



~a 



R M 
u~ 



Illustrative examples are presented in section 11. 



(70) 



11.0 Results and Conclusions 

An analysis for the prediction of the strength of fiber 
reinforced composite plates subjected to loads and thermal 
environment, described in this report, has resulted in a com- 
puter program for implementation of the analysis. A listing 
of the program, as well as directions for it’s use, is given 
in appendices B and C. 

Several analyses have been carried out in order to com- 
pare the results of the present analysis to other analytical 
results (ref. (4)), as well as experimental results. The fail- 
ure criterion used in Ref. (4) is the Tsai-Hill criterion. 
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Figures 6 through 8 present some preliminary results. Figure 
6 presents four sets of results for cross-ply laminates sub- 
jected to uniaxial loading: 

1) experimental data for E-glass/epoxy from Ref. (4), 

2) analytical results for E-glass/epoxy from Ref. (4), 

3) analytical results for E-glass/epoxy from the present 
analysis, and 

4) analytical results for T300/5208 graphite/ epoxy from 
the present analysis. 

The abscissa M, called the cross-ply ratio and defined by 

M = Et^/ Zt j for even j, and odd k (71) 
is a measure of the relative quantity of 90 and 0 degree 
plies. It is seen that for cross ply laminates under uniaxial 
loading, there is excellent agreement between experimental 
data and analytical results for the E-glass/epoxy laminates. 

All results are for a curing temperature of 270 degrees Fah- 
renfeit, and a loading environment at 70 degrees. In each case, 
the laminate did not fail at first ply failure. The first ply 
failure load is given in parentheses below the laminate fail- 
ure load. 

Figure 7 presents the same set of results for angle ply 
laminates as Figure 6 does for cross-ply laminates. In this 
case, the laminate has the layup with equal thickness 

plies. It can be seen that results from the present analysis 
for E glass/epoxy laminates are in good agreement with the 
experimental and analytical results of Ref. (4) for less 
than 30°, and <p greater than 60°. Tsai, Ref. (5), has a 
forthcoming report which shows that the Tsai-Wu criterion 
is in good agreement with experimental results for T300/5208 
graphite/epoxy angle ply laminates. 

Figure 8 shows the effect of temperature on laminate 
failure for t45° angle ply laminates under uniaxial loading 
for both E-glass/epoxy, and T300/5208 graphite epoxy laminates. 

In order to gain an understanding of the effect of 
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temperature on laminate strength, a series of computer 
analyses will be undertaken. The results of this computation- 
al effort will be presented in another report. 
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Appendix A 



Typical Properties for Two Lamina 

The table below presents some properties for two la- 
mina, E-glass/epoxy, and graphite/epoxy. The information 
is taken from references (3) and (4). 



Property 



Units 

SI (Engr'g.) 



E-glass/epoxy graphite/epoxy 



GPa 


(10 6 psi) 


53.8 


(7.8) 


181. 


(26.2) 


GPa 


(10 6 psi) 


17.9 


(2.6) 


10.3 


(1.49) 


dimensionless 


0.25 


(.25) 


0.28 


(.28) 


GPa 


(10 psi) 


8.6 


(1.25) 


7.17 


(1.04) 


dimensionless 


. 083 


(.083) 


0.16 


(0.16) 


1/C 


( 1/F) 


6 . 3y 


(3.5)y 


0. 02y 


( . 011) y 


1/C 


( 1/F) 


20.5 


(11. 4) y 


22. 5y 


(12. 5)y 


MPa 


(ksi) 


1034. 


(150) 


1500. 


(217.5) 


MPa 


(ksi) 


1034. 


(150) 


1500. 


(217.5) 


MPa 


(ksi) 


27.6 


(4.) 


40. 


(5.8) 


MPa 


(ksi) 


138. 


(20.) 


246. 


(35.7) 


MPa 


(ksi) 


41.4 


(6.) 


68 . 


(9.9) 


dimensionless 






0.7 


(0.7) 



v 



xy 

v 

Vx 



a 



a 



X 

X' 

Y 

Y' 

S 

+ v. 



* calculated by EQ.(12) 
t volume fraction of fiber 
y denotes micro (i.e.,10 ^) 
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Appendix B 



COMPUTER LISTING 
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oooooooooooooooooooooonoooooooooo 



FILE: STREN1 FORTRAN A 



NAVAL POSTGRADUATE SCHOOL 



C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 



AAAAAAAAAAAAtA A A A At At A A A A A A AAA A At At AA it AAA AAA* AAA A A A A A A A A A A At A** A* A 
* * 

* ANALYSIS OF THE STRENGTH OF A SYMMETRIC COMPOSITE * 

* PLATE DUE TO A THERMAL ENVIRONMENT * 

* * 
************************************************************ 



THIS IS THE MAIN PROGRAM WHICH READS IN THE PROPERTIES OF THE 
COMPOSITE, AS WELL AS THE INPLANE LOADS AND THERMAL DISTRIBU- 
TION THROUGH THE LAMINATE. IN ADDITION, THIS PROGRAM CALLS THE 
SUBROUTINE CYCLE, WHICH IN TURN CALLS SUBROUTINE FAILUR. 



STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 
STROO 

COMMON ANGLE ( 100 ) ,Z ( 101 ) , A ( 3 , 3 ) ,D( 3 , 3 ) , AT ( 3 , 3 ) , DT ( 3 , 3 ) , A I NV ( 3 , 3),WSTR00 
1K( 50),DINV( 3.31 ,PLYTHK( 100) , OB AR11 ( I 00) . QBAR12 { 100 ) ,QBAR22( 100) ,QBSTROO 
2AR16(100) ,QBAR26( 100 ) , QB AR66 ( 1 00 ) ♦ AL PHAl ( 100 ) , ALPHa2 ( 100) . AL PH A3 (1STR00 
300) ,CELTEM( 100) , STRES1 (100 ) , STRES2 ( 1 00 ) , STRES3 ( 100 ) , T EM ( 100 ) , STRESSTR 00 
4X( 100), STRESY(100),TAUXY( 1 00 ) , R 1< 100 ) ,R2 ( 100) , ZN1 , ZN2 ,ZN3 , FXX , FX ,FSTROO 
5YY,FY,FSS»FXY,QXX»QXY,QYY,QSS» ALFAX ,ALFAY,STFRTM,STRSL1(100 ) , STROO 

6STRSL2( 100) ,STRSL3(100) ,STRST1(100) ,STRST2(100) , STP.ST3 ( 100 ) , STROO 

7STRSTXC 100) ,STRSTY( 100 ) , TAUTXY ( 100 ) ,STRSLX( 100) ,STRSLY( 100) ,TAULXYSTROO 
8 ( 1 CO ) , R ( 70) ,RMAX , ZM1 ,ZM2 »ZM3 , STROO 

SNPLY,ICYCLE,IFAIL»NPL»IT( 100) STROO 

STROO 

AND 2 AFTER A NAME DENOTE THE QUANTITY IS N.R.T. THE PLATE STROO 

AXES (I.E..THE GLOBAL AXES OF THE SYSTEM). STROO 

X ANC Y DENOTE THE NATURAL COORDINATES OF A PLY. X IS ALONG THE STROO 

FIBER AND Y IS THE TRANSVERSE DIRECTION STROO 

A ANC C ARE THE INPLANE AND BENDING STIFFNESS MATRICES OF THE COMPO-STROO 
SITE PLATE W.R.T. LAMINATE AXES. THEY CHANGE AS PLIES FAIL. STROO 

A I NV AND D IN V ARE THE INVERSES OF THE A AND D MATRICES, RESP ECT IVELYSTROO 
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OF THERMAL EXPANSION W.R.T, 



THERMAL 



ALFAX, AND ALFAY ARE THE COEFFICIENTS 
THE NATURAL COORDINATES OF A PLY 
ALPHAl(I) ,ALPHA2(I ) ,ALPHA3(I ) ARE THE COEFFICIENTS OF 
EXPANSION OF THE I-TH PLY W.R.T. THE LAMINATE AXES. 

ANGLE ( I ) IS THE ANGLE BETWEEN THE FIBER DIRECTION ( I . E. ,X-AX IS ) OF 
THE I-TH PLY AND THE 1-AXIS OF THE LAMINATE COORDINATE SYSTEM. 

AT AND DT ARE THE A AND D MATRICES DIVIDED BY THE PLATE THICKNESS 
DELTEM(I) IS THE TEMPERATURE DIFFERENCE OF THE I-TH PLY, AND 
IS EQUAL TO TEM< I J-STFRTM. 

FXX,FX,FYY,FY,FSS ARE THE FAILURE PARAMETERS IN THE TSAI-WU FAILURE 
CRITERION AND ARE OBTAINED FROM X, XPRIME, Y, YPRIME, SS, AND 
FXYSTR. 

IT ( I ) EQUALS 0 FOR NONFAILED PLIES, AND EQUALS 1 FOR FAILED PLIES 
NPLY IS THE NUMBER OF PLIES ABOVE THE MIDPLANE OF A SYMMETRIC 
LAMINATE. THE TOTAL NUMBER OF PLIES ' ‘ 

PLIES DECREASES DURING THE ANALYSIS 
PLYTHK(I) IS THE PLY THICKNESS OF THE 
QXX,CXY,QYY» AND QSS ARE THE STIFFNESS 
THE LAMINA X AND Y AXES. 

QBAR11 (I ) , QBAR12 ( I ) ,QBAR22( I ),QBAR16( I ) , QB AR26 ( I ) , AND QBAR66 ( I ) 

ARE THE STIFFNESS COEFFICIENTS OF THE I-TH PLY W.R.T. THE PLATE 
(I.E., LAMINATE) COORDINATES. 

IS THE LOAD FACTOR (I.E. ,THE STRENGTH RATIO) FOR THE GIVEN 
TEMPERATURE DISTRIBUTION AND UNSCALED INPUT LOADING. R GREATER 
THAN 1. IS ASSOCIATED W. SAFE LOADING, AND R LESS THAN 1. MEANS 



IS 2*NPLY. THE NUMBER OF 
AS PLIES FAIL. 

I-TH PLY. 

COEFFICIENTS OF A PLY W.R.T. 



R 



STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 

STROO 
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FAILURE HAS GCCURRED. FAILURE OCCURS WHEN ABS <MIN<R<IIJ ) IS STR00560 

EQUAL TO UNITY. STR00570 

STRESl(I) »STRES2(I),STRES3(I) ARE THE STRESSES IN THE I-TH PLY W .R . T STR00580 
THE LAMINATE AXES. STR00590 

STRESX ( I ) » STRESYC I ) » AND TAUXY(I) ARE THE STRESSES IN THE I-TH PLY STR00600 
W.R.T. THE NATURAL (PLY) AXES, AFTER SCALING WITH THE LOAD FACTOR. STR 00610 
STRSL FOLLOWED 8Y A LAMINA OR LAMINATE '•SUBSCRIPT'' (I.E. STRSLX CR STR00620 
STRSL1) IS THE PLY STRESS DUE TO LOAD W.R.T. LAMINA OR LAMINATE STR00630 
COORDINATES BEFORE SCALING WITH THE LOAD FACTOR. STR00640 

STRST FOLLOWED BY A LAMINA OR LAMINATE "SUBSCRIPT" IS THE PLY STRESSSTR00650 



DUE TO TEMPERATURE W.R.T. LAMINA OR LAMINATE COORDINATES 

STFRTM IS THE STRESS FREE TEMPERATURE 

TEM(I) IS THE TEMPERATURE OF THE I-TH PLY. 

THICK IS THE PLATE THICKNESS 

XfXPRIMEtYtYPRIMEtSS, AND FXYSTR ARE THE FAILURE PARAMETERS DETER- 
MINED EXPERIMENTALLY. 

Z(I) AND Z(I+1) DENOTE THE LOCATIONS OF THE BOTTOM AND TOP 
SURFACES OF THE I-TH PLY, RESPECTIVELY. 

ZN1 , ZN2 , AND ZN3 ARE THE UNSCALED INPLANE LOADS 

ZNUX Y , ZNU YX , EX , EY , AND G ARE THE STIFFNESS COEFFICIENTS OF A PLY 
W.R.T. PLY AXES. 

READ (5,40 NPLY, ITYPE.ZNUXY, EX,EY,G» ALFAX.ALF AY, STFRTM 
ZNUYX=ZNUXY*EY/EX 
WRITE (6,50) NPL Y , I TYPE 
WRITE (6,60) 

WRITE (6,70) 

WRITE (6,71) 

WRITE (6,80) 

WRITE (6,90) 

WRITE (6,100) 

WR ITE (6, 105) 

NN=NPLY+1 
NPL-NPL Y 



( ZNUXY , ZNUYX, 
ALFAX, ALFAY 
STFRTM 



EX » EY ,G) 



IF I T YPE=0 ALL PLIES ARE THE SAME THICKNESS 
IF I TYPE® 1 PLIES HAVE DIFFERENT THICKNESSES 

ICYCLE=1 
Z( 1)=0. 

IF ( ITYPE.NE.O) GO TO 20 

READ (5,110) THICK 

PLT=THICK/NPLY 

DO 10 I =1 » NPLY 

IT ( I )=0 

PL YTFK( I ) = PLT 

Z( I+1)=Z(I )+PLT 

CONTINUE 

CONTINUE 

IF ( ITYPE.NE.l) GO TO 30 

READ (5,120) ( J,PLYTHK( I ) , 1=1, NPLY) 

DC 25 1=1, NPLY 

Z(I + 1)=Z(I)+PLYTHK( I) 

IT (I )=0 
CONTINUE 



STR00660 
STR00670 
STR00680 
STR00690 
STR00700 
STR00710 
STR 00720 
STR00730 
STR00740 
STR 00750 
STR00760 
STR00770 
STRC0730 
STR00790 
STR00800 
STR00810 
STR00820 
STR 00830 
STR00840 
STR00850 
STR 00860 
STR008 70 
STR00880 
STR00890 
STR00900 
STR009 10 
STR00920 
STR00930 
STR00940 
STR00950 
STR00960 
STR 00970 
STR00980 
STR00990 
STR 01000 
STR010 10 
STR01020 
STR01030 
STR01040 
STR 01050 
STR01060 
STR01070 
STR 01080 
STR01090 
STR01100 
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30 



C 

C 

c 



c 

c 

c 



c 

c 

c 



c 

c 



c 

c 

40 

50 

60 

70 

71 



RE 


AD (5 


, 130) 


( J f 


WR 


IT E ( 


6,65) 




WR 


ITE ( 


6 ,75 ) 


( J 9 


WR 


ITE ( 


6,85) 




WR 


ITE ( 


6,95) 


(Z< 


INPL 


T THE 


APPL 


I ED 


READ (5 


,140) 


ZN1 


READ( 5, 


140) 


ZMl , 


WR 


ITE(6 


,210) 




WR 


I TE ( 6 


,220) 


ZN1 


WR 


ITE< 6 


,230) 




WR 


ITE (6 


,240) 


ZMl 


WR 


I TE ( 6 


, 245) 





INPLANE LOADS 



LAMINA STIFFNESS COEFFICIENTS IN NATURAL (X,Y) COORDINATES 



CENCM=1 .0-ZNUXY*ZNUYX 

QXX=EX/DENOM 

CXY=(ZNUXY*EY)/DENQM 

CYY=EY/DENOM 

CSS=G 



WRITE (6,150) 

WRITE (6,160) QXX , GY Y, QXY , QSS 



CALCULATE THE COEFFICIENTS OF THE FAILURE CRITERION 



READ (5,170) X,XPRIME,Y,YPRIME,SS,FXYSTR 
WRITE(6,175) 

WRITE (6,180) X, XPRIME,Y,YPRIME,SS,FXYSTR 



FXX=l./( X*XPRIME) 

FX= (l./X)-( l./XPRIME) 
FYY=1 ./( Y*YPRIME ) 

FY = ( l./Y)-< l./YPRIME) 
FXY=FXYSTR*SQRT< FXX+FYY) 
FSS=l./( SS**2 ) 



WRITE (6,195) 

WRITE (6,190) 

WRITE ( 6,200) FXX»FX»FYY»FY»FSS»FXY 



CALL CYCLE 



RETURN 



FORMAT 

FORMAT 

FORMAT 

FORMAT 



( 2I3,7G10.4> 

( 1H0 , 2X , 2 8HNC . OF PLIES ABOVE M IDPLANE= , I 3 . 10X , 6HI 
( 1H0,2X,35HITYPE=0 FOR UNIFORM PLY THICKNESSES) 

( 2X» 39HITYPE= 1 FOR NON UNIFORM PLY THICKNESSES) 



FCRMAT( //1H0,2X ,' LAMINA ELASTIC PROPERTIES W.R.T. LAMINA 



STR01 
STROi 
STR01 
STROI 
STROI 
STROI 
STROI 
STROI 
STRCl 
STROI 
STRCl 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROi 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
STRCl 
STROI 
STROI 
STROI 
STROI 
STROI 
STROI 
TYPE=, 13) STROI 
STROI 
STROI 
COORD I NAT STROI 
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£ES : ' ) STR 01660 

FORMAT ( 1H0 ,2X , 5HNU XY= , G9. 3, 2X , 5HNUY X= , G 9. 3, 2X ,3HE X= , G9. 4 , 3X , 3HE Y=STRO 1670 
1 ,G9.4,3X,2HG=,G9.4) STR01680 

FORMAT ( IH0,2X,7HALPHAX=.G10.4,10X,7HALPHAY=,G10.4) STR01690 

FORMAT ( // 1H0 »2X » ’THE STRESS FREE TEMPERATURE IS',G10.4) SJR0L70O 

FORMAT! //1H0.2X, • LAMINATE CONSTRUCT I ON : • ) 

FORMAT !G10.4) 

( 13 ,G10.4 ) 

( I 3? 2G 10. 4 ) 

(3G10.4) 

( //1H0,2X ,46HSTIFFNESS COEFFICIENTS IN NATURAL 



FORMAT 

FORMAT 

FORMAT 

format 
1 ) 

FORMAT 



STR01710 
STR01720 
STR01730 
STR 01740 
ST RO 17 50 
COORDINATES: STR 01760 



STR01770 

( 1H0 *5X,4HQXX= *G10. 3,2X ,4HQYY=,G10. 3, 2X ,4HQXY= ,G 10. 3 ,2X ,4HQSTRC1780 
l SS= » G10 • 3 ) STR 01790 

FORMAT ( 6G10.4) STR01800 

F0RMAT(//1H0,2X, 'STRENGTHS W.R.T. LAMINA AXES:') ✓ STR01810 

FORMAT ( 1H0 , 4X » 2HX= ,G10 . 4, 3X , 7 HXPR IME= ,G 10 .4,3X ,2HY=, G10.4, 3X» 7HYPSTRC1820 
1RIME=,G10.4,3X,3HSS= »G10.4,3X» 7HFXYS TR= ,G10.4) STR01830 

FORMAT (1H0»7X*3HFXX,11X*2HFX,11X,3H FYY » 10X, 2HFY, 10X, 3HFSS, 1 OX , 3HF STR C 1840 



1XY) 

FCPMAT(//1H0,2X, 'TSAI-WU STRENGTH PARAMETERS W.R.T. LAMINA AXES:' 
FORMAT (4X,G10.4,3X ,G10.4,3X ,G10.4, 3 X, G1 0.4, 3X ,G 10 . 4, 3X, GIO .4 ) 
FORMAT! //1H0,2X,' INPLANE LOADS W.R.T. LAMINATE AXES:*) 

FORMAT! IH0,5X,*ZN1=' ,G10. 4 , 5X , ' ZN2=* ,G10.4, 5X, » ZN3= * ,G10. 4) 
FORMAT! //1H0.2X, • MOMENTS W.R.T. LAMINATE AXES:') 

FORMAT! IHOt 5X , • ZM 1= • ,G10. 4 , 5X , • ZM2 = ' ,G10 .4 ,5 X , • ZM3 = • ,G10. 4 ) 
FORMAT! 1H0.2X,'*** MULTIPLICATION OF THE A80VE LOADS AND MOMENTS 
1Y THE STRENGTH RATIO AT LAMINATE FAILURE GIVES THE FAILURE LOADS 
2ANC MOMENTS * ) 

FORMAT ( 1H0 » 9X » 3HPLY ,10X ,5HANGLE,11X , 13HPLY TH ICKN ESS,7X , 3HTEM ) 

! 9X, I3,10X,G10.4,9X,G10.4,6X,G10.4) 

! IH0.2X.10HZ- LOCATION) 

(5X.G10.4) 



FORMAT 

FORMAT 

FORMAT 

END 



SUBROUTINE CYCLE 

SUBROUTINE CYCLE CONSTRUCTS THE A AND D MATRICES, AND THEN SOLVES 
FOR THE STRESSES IN EACH PLY. AFTER THE STRESSES ARE CALCULATED, 
SUBROUTINE FAILUR IS CALLED TO DETERMINE WHICH, IF ANY, PLIES 
HAVE FAILED. SUBROUTINE FAILUR IN TURN CALLS SUBROUTINE CYCLE 
UNTIL ALL PLIES WHICH WILL FAIL DUE TO THE LOAD AND THERMAL 
ENVIRONMENT HAVE FAILED. 

COMMON ANGLE ( 100 ),Z! 101) .A (3,3 > ,D( 3. 3 ) , AT! 3,3 ) ,0T ! 3 ,3 ) . A I NV < 3 , 3 ) 



STR 01850 
) STR01860 
STR01870 
STR01880 
STR 01 890 
STR 01900 
STR01910 
B STR 01920 
STR01930 
STR01940 
STR 019 50 
STR 019 60 
ST RO 19 70 
STR 01980 
ST RO 1990 
STR02000 
STR02010 
STR02020 
STR 02030 
STR 02040 
STR02050 
STR02060 
STR02070 
STR02080 
STR 02090 
STR02100 
WSTR 02110 



turriLIM Ai!^LC IAUU I i'J 4 I |A U ) J I f U \ J J O I f A I \ Jf u i f LM l jMrtillV ! i WO 

1K( 50) ,DINV( 3,3) .PLYTHK(IOO) ,QBARll(100) , QBAR12(100) ,Q8AR22l 100 ) , Q8STR02 120 
2AR16(10C) ,Q BAR 2 6 ( 100J.QBAR66! 100) .ALPHA! ( 100) ,ALPHA2(100) , ALPHA3 ! 1STR02130 
300) ,CELTEM(100) , STRESldOO ) , STRESS ( TOO ) , STRES3 ( 1 00 ) , TEM ( 100 ) , STRESSTR02140 
4X( ICO) , STRE SY( 100 ) , TAUXY ( 1 00 I , Rl( 100 > ,R2 < 100 ) , ZNl , ZN2 »ZN3 , FXX , FX , FSTR021 50 
5YY,FY.FSS.FXY,QXX,3XY,QYY,QSS, ALFAX, ALFA Y,STFRTM. STRSL1! 100) , STR02160 

6STRSL2 ( 1 00) ,STRSL3(100) ,STRST1 ( 100 ) , STRS T2 ( 100 ) , STRST3! 100 ) , STR 02 170 

7STRSTX! 100) ,STRSTY( 100) ,TAUTXY!100) ,STRSLX(100) ,STRSLY(100) »T AULXYSTR02180 



8(100), R(70) ,RMAX • ZM1 , ZM2 , ZM3 » 
9NPLY,ICYCLE,IFAIL,NPL,IT( 100) 



STR 02 1 90 
STRC2200 
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C 

C 

c 



c 

c 

c 



c 

c 



WRITE (6,50) 
NN = NPL + 1 



ICYCLE 



ZERO OL'T THE A ANO D MATRICES, AND THE THERMAL LOAD VECTOR 



A(l, 1) = 0. 
A< 1*21=0. 
A ( I, 2 )=0 . 



A ( 2 « 2 ) =0 • 
At 2, 31 = 0. 



A ( 3 , 3 )=0 • 

0 ( 1 , 1 )= 0 . 



D( 1, 2) = 0. 

C C 1 , 31=0 . 



D(2,2)=0. 

C( 2, 3 ) = 0 • 

D(3,3)=0. 

ZNT1=0. 

ZNT2=0. 

ZNT3=0. 



LAMINA STIFFNESS COEFFICIENTS W.R.T. LAMINATE COORDINATES 



WRITE 

WRITE 



( 6 , 100 ) 

( 6 , 110 ) 



CO 20 IPLY=1,NPLY 
IF ( IT ( IPLY ) . EQ . 1 ) GO TO 20 
DELTEM( IPLY ) = TEM ( IPLY )-STF RT M 
ANGLE1=ANGLE( IPLYJ/57.2957 8 
S= S I N ( ANGLE1 ) 

S2 = S**2 
S3=S*S2 
S4=S2**2 
C = CCS (ANGLE 1 ) 

C2=C**2 
C3=C2*C 
C4=C 2**2 
CSSQ=C2*S2 

QBAR1K IPLY) = CXX*C4+2 . *( QXY+2. *QSS ) *CSSQ+QYY *S4 
QBAR12( IPLY )= ( QXX+QYY-4. *Q SS) *CSSQ+Q XY*< S4+C4J 
QB AR22( IPLY)=CXX*S4+2.*(QXY+2.*QSS)*CSSQ+QYY*C4 
QBAR16( I PLY ) = ( QXX-QXY-2. *QSS ) *S*C3+( QXY-QYY+2 . *QSS ) *S3*C 
CBAR26( IPLY1=(QXX-QXY-2.*QSS )*S3*C+< QXY-QYY+2. *QSS) *S*C3 
QBAR66 ( IPLY ) = ( QXX+QYY— 2 . *QXY-2 • *QSS ) *CSSQ+QSS* ( S4+C4 ) 



ALPHA1 ( IPLY J=C2*ALFAX+S2*ALFAY 
ALPHA2 ( IPLY1=S2*ALFAX+C2*ALFAY 
ALPHA3( IPLY|=(ALFAY-ALFAX) *C*S 



WRITE (6,120) IPLY, QBAR1K IPLY), QBAR22( IPLY), QBAR12( IPLY) 
IPLY) »QBAR16(IPLY) ,QB AR26( I PLY ) 



CONSTRUCT THE THERMAL LOAD VECTOR, ZNT 



STR022 
STR022 
STR022 
STR022 
STR022 
STR022 
STR022 
STR022 
STR022 
STR023 
STR 023 
STR023 
STRC23 
STR023 
STR023 
STR023 
STR023 
STR 023 
STR 023 
STR024 
STR 024 
STR 024 
STR024 
STR024 
STR 024 
STR024 
STR024 
STR 024 
STR024 
STR025 
5TR025 
STR025 
STR 025 
STR025 
STR025 
STR 025 
STR025 
STR025 
STR025 
STR026 
STR026 
STR026 
STR 02 6 
STR 026 
STR026 
STR 026 
STR 026 
STR026 
STR026 
STR027 
STR 027 
,QBAR66( ISTR027 
STR027 
STR 027 
STR 027 
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P = P LYTHK (IPLY) 

DTM=CELTEM( IPLY) 

A1=ALPHA1(IPLY) 

A2=ALPHA2(IPLY) 

A2=ALPHA2( IPLY) 

ZNT 1=ZNT1 + P*DTM*(QBAR11(IPLY)*A1+QB-AR12( IPLY) * A2+QBAR16( IPLY)*A3 ) 
ZNT2 = ZNT2«-P*DTM*(C3AR12( IPLY )*A1+QBAR22 < IPLY )*A2+QBAR26( IPLY)*A3 ) 
ZNT3=ZNT3+P*DTM*( QBAR16( IPLY)*AH-QBAR26( IPLY)*A2*QBAR66 ( 



IPLY )*A2+QBAR66 ( IPLY )* A3) 



CONSTPUCTICN OF LAMINATE A AND D MATRICES FROM LAMINA MATRICES 

Z1=Z ( IPLY+l ) 

Z2=Z ( IPLY ) 

THICK=2 • ♦ ( Z 1 — Z 2 ) 

RK=2 • * ( (Zl**3 )-( Z2**3) )/3. 



A ( 1 , 1 )= 
A( 1,2) = 
A( l, 2 ) s 
A ( 2 , 1 ) = 
A( 2,2) = 
A( 2,2 ) = 
A (3 ,1 ) = 
A( 3,2) = 
A ( 3 » 3 ) = 

DC 1,1) = 
C( 1,2) = 
D C l ,3) = 
D( 2,1)= 
0(2,2) 
0 ( 2 , 2 ) = 
C(3, l) = 
0(3,2) = 
0(3,2) = 



A(1,1)+THICK*QBAR11(IPLY) 
=A ( 1 ,2 ) +THICK*CBAR12 (IPLY) 
= A ( l » 3 ) + TH I CK*QBAR 16(IPLY) 
= A ( 1 , 2 ) 

=A(2,2)+THICK*QBAR22 (IPLY) 
=A(2,3)+THICK*QBAR26(IPLY) 
=A ( 1 , 3 ) 

=A ( 2 , 3 ) 

A(3,3)* , THICK*QBAR66(IPLY) 



= D( l, 
= 0 ( 1 , 
= 0 ( 1 , 
= D ( 1 » 
= 0 ( 2 , 
= 0 ( 2 , 
= 0 ( 1 , 
= 0 ( 2 , 
= D ( 3 , 



1 ) +RK*QBARll ( IPLY) 
2 >+RK*C3AR12( IPLY ) 
3)+RK*C8ARlo(IPLY) 

2 ) 

2 >+RK*CBAR22( IPLY ) 
3) +RK*QBAR26( IPLY) 
3) 

3 ) 

3)+RK*QBAR66( IPLY) 



CONTINUE 



ZNT 1=2. *ZNT 1 
ZNT2=2. *ZNT2 
ZNT2=2.*ZNT3 

CALCULATE THE TOTAL INPLANE FORCES. T0TN1= ZN1+ZNT1 

T0TM = ZM* , ZNT1 
TCTN2=ZN2+ZNT2 
T0TN2=ZN3+ZNT3 

UR ITE( 6, 125) 



WRITE 

WRITE 

WRITE 



(6 



ZN1 ,ZNT1,T0TN1 

ZN2,ZNT2,T0TN2 

(6,150) ZN3,ZNT3,T0TN3 



(6 



,130) 

140) 



WRITE (6,155) 



STR02760 
STRC2770 
STR 02780 
STR02790 
STR02800 
STR 028 10 
STR 02820 
STR02830 
STR02840 
STR 02 850 
STR 02860 
STR02870 
STR 02 8 80 
STR 02890 
S-TR02900 
STR029 10 
STR02920 
STR02930 
S TR 02940 
STR02950 
STR 02 9 60 
STR 02970 
STR02980 
STR 02990 
STR 03000 
STR03010 
STR03020 
STR03030 
STR 03040 
STR 03050 
STR03060 
STR 03070 
STR 03080 
STR03090 
STR03100 
STR03110 
STR03120 
STR 03130 
STR03140 
ST^03150 
STR03160 
STR03170 
STR03180 
STR03190 
STR03200 
STR 032 10 
STR03220 
STR03230 
STR 03240 
STR03250 
STR03260 
STR 03270 
STR03280 
STR 03 2 90 
STR 03300 



39 



noon o oonoooo o ooo 



FILE 



STREN1 FORTRAN A NAVAL POSTGRADUATE SCHOOL 



25 

C 



WRITE (6,160) 

DO 25 I = 1 » N PL Y 

IF ( IT ( I ) . EQ . 1 ) GO TO 25 

WRITE (6,170) I,ALPHA1(I) ,ALPHA2(I) , ALP HA 3(1) 
CONTINUE 



WRITE (6,180) 
WRITE 



(6,135) A( 1,1 ) ,A(1 ,2 ) ,A( 1,3) ,D (1, 1 ),D(1,2) ,D( 1,3) 
WRITE (6,190) A(2,l) ,A(2, 2), A( 2,3) ,0(2,1) ,0(2,2) ,0(2,3) 



WRITE (6,190) A(3,l > ,A(3,2 ),A(3,3) ,0(3, 1 ),0(3,2) ,0(3,3) 
THE A AND C MATRICES PER UNIT THICKNESS OF LAMINATE 



H=2. *Z(NN) 

WRITE (6,200) H 
AT ( 1,1) =A( 1, 1)/H 
AT ( 1 , 2 ) =A( 1 , 2 ) /H 
AT ( 1 ,3 ) ®A ( 1 ,3)/H 
AT ( 2 • 1 ) =AT( 1,2) 
AT (2»2)=A(2,2)/H 
AT (2 ,3) =A(2 ,3)/H 
AT( 3 , 1) =AT( 1,3) 
AT (3,2) a AT ( 2 ,3) 
AT(3,3)=A(3,3)/H 



OT ( 1 ,1) a 0( 1 , 1 )/H 
0T( 1,2) =C( 1,2)/H 
CT ( 1 ,3) -0( 1 , 3 ) /H 
0T(2,1)=DT( 1,2) 

CT ( 2 » 2 ) -D( 2,2)/H 
CT(2,3)=0(2,3)/H 
CT ( 3 ,1 ) =0 T( 1 ,3 ) 

DT ( 3,2)=0T(2,3) 

OT (3 ,3) =C( 3 ,3 )/H 
WRITE (6,210) 

WRITE (6,185) AT ( 1 , 1), AT (1,2) ,AT(1,3) ,DT(1,1) , OT ( 1 , 2 ) ,0T(1,3) 
WRITE (6,190) AT ( 2 , 1 ),AT (2, 2), AT (2, 3 ),0T (2,1),0T(2,2),0T(2,3) 
WRITE (6,190) AT ( 3, 1 ) , AT ( 3 ,2 ) » AT(3 ,3 ) ,0T (3,1) ,DT(3,2),0T(3, 3) 



INVERSION CF THE A AND 0 MATRICES 



CALL LINV2F ( A, 3 , 3, A INV, 5 , WK , I ER) 
CALL LINV2F (0,3 ,3, DINV,5, WK, I ER) 



WRITE | 6 ,220 ) 



185) AINV( 1 ,1) ,AINV (1 ,2),AI NV( 1 ,3 ) ,0 INV (1,1),0INV( 1,2) 



WRITE 
1NV( 1,3) 

WRITE (6,190) AINV(2,1),AINV(2,2),AINV(2,3),0INV(2,1),0INV(2,2) 
INV (2,3) 

WR ITE (6,190) AINV( 3,1) ,AINV( 3,2) ,AI NV(3 ,3) ,01NV(3,1) ,DINV(3,2 ) 
INV (3 ,3) 



CONSTRUCT THE EXPRESSION FOR PLY STRESSES DUE TO TEMPERATURE 
1. A INV*ZNT 



STR033 
STR 033 
STR033 
STR033 
STR 033 
STR033 
STR033 
STR 033 
STR033 
STR 034 
STR 034 
STR034 
STR034 
STR034 
STR034 
STR 034 
STR034 
STR034 
STR 034 
STR035 
STR035 
STR035 
STR035 
STR 035 
STR035 
STR035 
STR 035 
STR035 
STR035 
STR036 
STR036 
STR 036 
STR036 
STR036 
STR 036 
STR036 
STR 036 
STR036 
STR036 
STR037 
STR037 
STR037 
STR037 
STR037 
STR037 
, 0 ISTP. C37 
STR037 
,01 STR037 
STR 037 
, DISTR038 
STR038 
STR038 
STR038 
STR 03 8 
STR038 
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WRI TE ( 6 » 225 ) 

WR I TE < 6, 230) 

DO 30 I = 1 »NPLY 

IF ( I T ( I ) .EQ.l ) GO TO 30 

AINT 1=A INV( I* 1) *ZNT 1 + AINV( 1, 2) *ZNT2+ AINV (1 ,3)*ZNT3 
AINT2=AINV(2,1)*ZNT1+AINV( 2t 2 ) *ZNT 2 + A INV ( 2 , 3 1 *ZNT3 
AINT3=AINV( 3,1)*ZNT1+AINV(3,2)*ZNT2+AINV (3,3)»ZNT3 

2. ALPHA*DELTEM 

ALTM1=ALPHA1 (I)*DELTEM( I) 

ALTM2=ALPHA2(I)*DELTEM( I) 

ALTM2=ALPHA 3 ( I ) *DELTEM( I ) 

3. THERMAL STRESS= GBAR*( AI NV*ZNT-ALPHA*DELTEM ) 

Al= A INT 1-ALTM l 
A2=AINT2-ALTM2 
A3=A INT3-ALTM3 

STRSTK I )=QBAR11 ( I)*A1+QBAR12 ( I )*A2*QBAR 16( I ) *A3 
STRST2( I ) =QBAR12 ( I) * A1+QBAR22 ( I j*A2 +QBAR26 ( I ) * A3 
STRST3C I )=QBAR16( I ) *A1+QBAR26( I )*A2*QBAR66( I ) *A3 
WRITE (6,240) I , STRSTK I ) »STRST2( I ) , ST R ST 3(1 ) , DELTEM( I ) 
CONTINUE 

CALCULATE THE PLY STRESSES DUE TO LOADS? STRSL=Q*AINV*ZN 

DO 35 1=1 ,NPLY 

IF ( IT( I ) .EQ.l) GO TO 35 

AINL1=AINV(1,1)*ZN1+AINV(1 ,2)*ZN2+AI NV( 1,3)*ZN3 
AINL2=AINV( 2,1) *ZNl+ AINV (2,2) *ZN2+AI NV(2,3J*ZN3 
AINL3=AINV(3,1)*ZN1 + AINV(3,2)*ZN2+AI NV( 3,3)*ZN3 

STRSLK I )=QBAR1 1 ( I ) * AINL 1+QBAR 12( I ) * A INL2+QBAR16 ( I ) *AI NL3 
ST RSL2 ( I ) = QBAR12 ( I) *A INL1+QBAR22( I ) *A INL2+QBAR26( I ) *AINL3 
STRSL3( I ) =QBAR16 ( I) * AINL1+QB AR26 ( I ) * A INL 2+CBAR66 ( I ) *A INL3 

CONTINUE 



CALCULATE THE STRESSES DUE TO MOMENTS 

IF( ZM 1. EC. 0. .AND . ZM2 . EQ.O . .AND.ZM3 .EQ .0 . ) GO TO 37 

DO 36 I =l,NPLY 

IF ( IT ( I ) • EQ. 1 ) GO TO 36 

ZZ=(Z(I+1)+Z( I) )/2. 

ZKAP 1=D INV( l,l)*ZMl+OINV(l ,2 ) *ZM2+DI NV( 1 ,3 ) *ZM3 
ZKAP2=DINV(2, 1 ) *ZM1 +C INV( 2 ,2 ) *ZM2+DI NV( 2,31 *ZM3 
ZKAP3=DINV(3,1) * ZM1 + C INV (3 ,2 )* ZM2+D I NV ( 3 , 3 ) *ZM3 

DINM1=ZZ*ZKAP1 
DI NM2=ZZ*ZKAP2 
DINM2=ZZ*ZKAP3 



STR03860 
STR03870 
STR03880 
STR03890 
STR03900 
STR039 10 
STR03920 
STR 03930 
STR03940 
STR03950 
STR 03960 
STR03970 
STR039 80 
STR 03990 
STR04000 
STR04010 
STR04020 
STR04030 
STR04040 
STR04050 
STR04060 
STR 04070 
STR04080 
STR04090 
STR04100 
STR041 10 
STR04120 
STR 04 130 
STR04140 
STR 04150 
STR04160 
STR 041 70 
STR04180 
STR04190 
STR04200 
STR04210 
STR04220 
STR 04230 
STR04240 
STRC4250 
STR04260 
STR04270 
STR 04280 
STR 04290 
STRC4300 
STR 043 10 
STR04320 
STR 04330 
STR04340 
STR04350 
STR04360 
STR04370 
STR043 80 
STR 04390 
STR04400 
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C 

C 

c 

c 



36 

37 



38 

C 

C 

c 

c 



c 

c 

c 

c 



c 

c 

c 



ADD TC STRESSES DUE TO INPLANE LOADS 



STRSL1 ( I ) =STRSL1 ( II+CBAR11 ( I ) *D INM1+QBAR 12 ( II *D INM2+QBAR 16 ( I) 
63 

STRSL2( I )=STRSL2 I I) +QBAR12 ( I ) *D INM 1+QBAR 221 I) *D INM2+QBAR26I I) 
63 

STR$L3( I > = STRSL3( D+QBAR 16 ( I >*DINM1+QBAR26( I ) *DI NM2+CBAR66 ( I) 
63 



CONTINUE 
CONTINUE 
IkR ITE(6 t 270 ) 



DC 38 I- ltNPLY 
IF ( I T ( I I .EQ.l) GO TO 38 



WRITEI6»280) 
CONTINUE 



I» STRSLKII f STRSL2I I ) t STRSL3II) 



FAILURE ANALYSIS 

9 1. CONVERT PLY STRESSES TO NATURAL (LAMINAI AXES 
DO 40 1= ltNPLY 
IF ( IT ( I I.EQ.l) GO TO 40 



ANGLE 1=ANGLE(I)/ 57. 29578 
S=SI MANGLE 1 ) 

S2=S=»*2 
C=COSI ANGLE1 ) 

C2=C**2 

SC=S*C 

STRSTXI I >=SIRSTli II*C2+STRST2( I )*S2+STRST3 (I)*(2.*SC) 
STRSTYI I J = STRST1< I) *S2+ STR ST2 (I) *C2- STRS T3 I I )*<2.*SC) 
TAUTXY ( I ) = SC* ( STRST2 ( I J-STRSTl ( I H+STRST3I I )*( C2-S2 ) 



STRSLXI I )=STRSL 1 I I> *C2+STRSL2( I ) *S2* S TRS L3 ( I ) * ( 2. * SC ) 



STRSLYI I J = STRSL1 ( I ) *S2+S TRSL2 ( I )*C2-STRSL3 ( I ) * { 2 .* SC I 
TAULXYI I» = SC*(STRSL2(I>-STRSLl(in+. r ' 



■STRSL3II >*( C2-S2 ) 



2. CALCULATE THE COEFFICIENTS OF THE QUADRAT IC:A*(R(I)**2) 
♦6*R(I)+C = 0. 



AA=FXX*( STRSLXI I )**2 )+2 . *F XY* ( STRSLX I I)*STRSLY II) >+FYY*(STRSL 
1*2I+FS$*(TAULXY( I)**2) 

B»2.*(FXX* STRSLXI I)*STRSTX(I ) +FXY*ST RSL X 1 1 )* STRSTYI I ) + 
1FXY*STRSLY(I)*STRSTX(I)+FYY*STRSLYII I *ST RSTY 1 1 ) +F S S *TAULXY I I) 
2*TAUTXYI I) ) +FX*STRSLX( I } -*-FY*STRSLY 1 1 ) 

C=FXX*( STRSTXII )**2 f +2. *FXY*STRSTX 1 1 ) *ST RSTY 1 1 } +FYY*< STRSTY ( I 
1+FSSMTAUTXYI I ) **2) + FX*$TR ST XI I J+FY* STRSTYI I | — 1 • 



3. SOLVE THE FAILURE QUADRATIC FOR THE STRENGTH RATIOS* R(I) 



CC=SCRT( (8**2 I -4 • *AA*C ) 
R 1 1 I )=(-B+CC>/I2.*AA) 
R2(n=(-e-CCI/(2.*AA) 



STR044 
STR 044 
STR044 
*D I NMSTR 044 
STR044 
*DI NMSTR044 
STR 044 
*DI NMSTR044 
STR044 
STR 045 
STR045 
STR045 
STR045 
STR045 
STR 045 
STR045 
STR 045 
STR 045 
STR045 
STR046 
STR046 
STR046 
STR 046 
STR046 
STR 046 
STR 046 
STR046 
STR 046 
STR046 
STR047 
STR 047 
STR047 
STR 047 
STR 047 
STR047 
STR047 
STR047 
STR047 
STR 047 
STR 048 
STR 048 
Y I I )*S7R048 
STR048 
STR048 
STR048 
STR048 
I **2 ) STR 048 
STR048 
STR 048 
STR 049 
STR049 
STR C49 
STR049 
STR049 
STR 049 
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GO TO 41 
I , STRSTXC I ) 



, STR STY ( I ) ,TAUTXY ( I ) ,R1II),R2(II 



CONTINUE STR04960 

STR04970 
STR 04980 
STR049 90 
STR 05000 
STR 050 10 
STR05020 
STR 05030 
STR05040 
STR05050 
STR 05060 
STR05070 
STR 050 80 
STR 05090 
STR 05 100 
STR05110 
STR 05 1 20 
STR05130 
STR 051 40 
AXES : STR05 150 
STR05160 

I 1H0,3X,3HPLY,3X,6HQBAR11,5X,6HQBAR22,5X,6HQBAR12,5X,6HQBARSTR05170 



WRITE ( 6 , 250 ) 
00 4 *-i-=-r,NPLY 
IF ( IT ( I ) • EQ .1 ) 
WRITE (6,260) 
CONTINUE 

WR ITE(o, 290) 

CC 42 I=1,NPLY 
IF( IT(I ).EQ.l) 
WR ITE( 6 ,300 ) 
CONTINUE 

CALL FAILUR 



GO TO 42 

I , STRSLXt I ) , STRSLYI I ) ,TAULXY(I) 



RETURN 

FORMAT 
FORMAT 
1 * ) 

FORMAT 



(1H1 ,1X,7HICYCLE=,I3) 

( //1H0,2X, 'PLY STIFFNESS COEFFICIENTS W.R.T. LAMINATE 



( 6X,4HZN3=,G10.4,10X,5HZNT3=,G10.4,10X,6HTOTN3=,G10.4) 

(//1H0,2X, ‘THERMAL COEFFICIENTS OF EXPANSION W.R.T. LAMINATSTR05240 



) 



166 , 5X,6HQ8AR 16, 5X,6HQBAR26) 

FORMAT (2X,I3,6( IX, G10 . 3 ) ) 

FQRMAT( // 1H0, 2X, 'APPLIED, THERMAL, AND TOTAL INPLANE FORCES:') 
FORMAT ( 1H0, 5X» 4HZN 1= ,G 1 0. 4, 10X, 5HZN Tl= ,G10.4 , 10X , 6HTCTN1* , G 10 .4 ) 
FORMAT (6X,4HZN2=,G10.4,10X,5HZNT2=, G10.4, 10X, 6HT0TN2=,G10. 4) 

FORMAT -* “* ‘ 

FORMAT 
IE AXES: 

FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT( 



STR05180 
STR051 90 
STR 05200 
STR05210 
STR05220 
STR05230 



STR05250 

( 1H0,3X,3HPLY,11X,6HALPHA1,14X,6HALPHA2 , 14X , 6 FALPHA6 ) STR05260 

(3X,I3,10X,G10.4,10X,G10.4»10X,G10.4) STR05270 

( // 1H0 » 12X » 8HA— MATR IX ,3 5X , 8 HD-MAT RIX ) STR 05280 

I 1H0, IX, Gil. 4, IX, Gil. 4, IX, Gil. 4, 8X, 311. 4, IX, Gil. 4, IX ,G11 ,4)STR05290 
(2X,G11.4, IX, Gil. 4, 1X,G11.4,8X,G1 1.4, IX, G1 1.4, IX, Gil. 4) STR05300 

( //1H0,5X,16HPLATE THICKNESS* ,G10 .4) STR05310 

</lH0,9X, 18HA-MATR IX/ THICKNESS, 25 X, 18HD-MATR I X/THI CKNESS) STR05320 
( 1H0 , 9X , 16HA-I NVERS E MATRIX, 30X, 16HD-INVERSE MATRIX) STR05330 

//1H0,2X, 'STRESSES IN LAMINA AND LAMINATE COORDINATES BEF0RSTR05340 
&E MULTIPLICATION WITH THE LOAD FACTOR:') STR05350 

FORMAT ( 1H0,3X,3HPLY,6X,6HSTRST1,12X,6HSTRST2, 14X , 6HSTRST3 , 9X , 6HDESTR05360 
1LTEM ) STR05370 

FORMAT (4X,I3 ,7X ,G10 .4,8X, G10.4,9X,G10.4,5X,G10.4) STR05380 

FORMAT ( 1H0 ,3X,3HPLY, 6X» 6H STRSTX,9X , 6HST RSTY , 1 IX » 6HT AUTXY , 1 2X , 2HR 1STR 05390 



, 5X.G10.4) 
6 HSTRSL3 I 



1, 12X.2HR21 

FORMAT ( 4X, 13 ,5X , G1 0 .4 ,5X , G1 0 . 4 , 5X , G 10.4 , 5X,G1 0 .4 
FORMAT( 1H0, 3X , 3HPL Y , 6X ,6HS TR SL 1 , 12X , 6HSTRSL2 ,14X, 

FORMAT ( 4X,I3,7X,G10.4,8X,G10.4,9X,G10.4) 

FORMAT( 1H0, 3X ,3HPLY , 6X ,6HSTRSLX ,9X , 6 HST RSLY , 1 IX , 6HT AULXY ) 
FORMAT ( 4X,I3, 5X,G 10.4, 5X,G 10.4, 5X,G1 0.4) 

END 

SUBROUTINE FAILUR 



STR05400 
STR 054 10 
STR05420 
STR 05430 
STR 05440 
STR05450 
STR05460 
STR 05470 
STR05480 
STR05490 
STR 05 5 00 
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C 

C 

C 

C 

C 

C 

C 

c 

c 

c 



10 

c 



c 

c 

c 



20 

C 



30 

c 

c 

c 



THIS SUBROUTINE DETERMINES WHICH PLIES HAVE FAILED BY USING 
THE TSAI-WU FAILURE CRITERION. THE PLY WHICH HAS FAILED IS 
ELIMINATED AND THE ANALYSIS IS RETURNED TO SUBROUTINE CYCLE, 
WHICH RECALCULATES THE STRESSES FOR THE MODIFIED LAMINATE AND 
RETURNS TO SUBROUTINE FAILUR AGAIN. THIS PROCESS IS REPEATED 
UNTIL ALL PLIES WHICH ARE TO FAIL HAVE FAILED. 



STR05S 
STR051 
STR055 
STR055 
STR05S 
STR055 
STR05E 
SIR 05f 
STR05 5 
STR054 

COMMON ANGLE! 100 ) . Z ( 101 ) ,A (3 ,3 ) ,D< 3 , 3 ) , AT ( 3 , 3 ) , DT ( 3 , 3 I , A INV ( 3, 3 J , WSTR 056 
IK(50),DINV(3,3) »PLYTHK! 100 ) , QB AR 11 ( 1 00) ,QBAR12( 100) ,CBAR22( 100) ,QBSTR056 
2AR16(10C) ,QBAR26 (100 ) ,QBAR66 ( 100 ) , AL PHA 1 (100 ) , ALPHA2 ( 100) .ALPHA 3 ( 1 STR 05* 
300) , DEL TEM{ 100) , STRE SI < 1 00 ) , STRES2 ( 1 00) ,STRES3 (100 ) » TEM( 100 ), STRESS TR05£ 
AX ( 100),STRESY(100), TAUXY ( 100 ) , R 1( 100 ) ,R2 ( 100) ,ZN1,ZN2,ZN3 ,F XX , FX ,FSTR05* 
5YY»FY,FSS,FXY »QXX ,QXY , QYY » QSS » ALFAX , ALFAY, STFRTM , STRSL1 ( 100 ) , STR05* 

6STRSL2( 100) ,STRSL3( 100) ,STRST1(100) , STRST2(100) ,STRST3(100) , STR05* 

7 STRSTX ( 100) ,STRSTY( 100),TAUTXY( 100 ) ,STRSLX( 100) ,STRSLY( 100) , TAULXYSTR05* 



THIS SUBROUTINE CHECKS FOR PLY FAILURE 



8( 100) ,R( 70) , RMAX , ZM1 ,ZM2.ZM3 , 
SNPLY, ICYCLE, I FAI L,NPL, IT l 100) 



R( ICYCLES = l.E+20 



DO 10 1=1 » NPLY 

IF ( I T( I ).EQ.l) GO TO 10 

A8SR1 = A£S (R 1 ( I ) ) 

ABSR2=ABS(R2(I) ) 

ABSR=AMIN1( ABSR1, ABSR2) 

IF <ABSR.LT.R( ICYCLE) ) IFAIL=I 

IF (ABSR.LT.R(ICYCLE) ) R( I C YCLE )=ABSR 

CONTINUE 



WRITE (6, 110) ICYCLE, I FAIL,R( ICYCLE) 



CALCULATE THE STRESSES AT IMPENDING PLY FAILURE 



CO 20 1=1, NPLY 

IF ( IT( I J.EQ.l) GO TO 20 

STRESX! I)=R( ICYCLE)*STRSLX(I )+STRSTX (I) 

ST RESY ( I)=R( ICYCLE) *STRSLY (I >+STRSTY( I) 

TAUXY (I )=R( ICYCLE )*TAULXY( I)+TAUTXY( I) 

CONTINUE 



WRITE(6, 100) 

WR IT E ( 6 , 120) 

DC 3C 1=1, NPLY 
IF ( IT( I J.EQ.l) 
WRITE 16, 125) ( 

CONTINUE 



GO TO 30 

I, STRESX( I ) , ST RESY (I ) ,TAUXY( I) ) 



CALCULATE THE STRENGTH RATIO FOR LAMINATE FAILURE 



IF ( ICYCLE. EC. 1) RMAX=R(1) 
RMAX =AMAX1(R( ICYCLE), RMAX) 



STR 05* 
STR05T 
STR051 
STR 051 
STR051 
STR051 
STR 051 
STR051 
STR 051 
STR051 
STR051 
STR 056 
STR056 
STR056 
STR 056 
STR056 
STR056 
STR 056 
STR056 
STR 056 
STR056 
STR051 
STRC51 
STR 051 
STR 051 
STR 051 
STR051 
STR051 
STR 051 
STR051 
STR051 
STR06C 
STR06C 
STR06C 
STR06C 
STR06C 
STR06C 
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ECK FOR LAMINATE FAILURE 


STR06060 




STR06070 


NPL=NPL-1 


STR 06080 


IF(NPL.NE.O) GO TO 50 


STR06090 


IaR ITE ( 6 * 130 ) RMAX 


STR061C0 


IFINPL.EC.O) STOP l 


STR 061 10 


ICYCLE = ICYCL E+l 


STR06120 


IT ( I FAI L ) = 1 


STR 061 30 




STR 06 140 


CALL CYCLE 


STR06150 




STR 06 160 


RETURN 


STR06170 




STR06180 


F0RMATC//1H0.2X,* STRESSES AT 


IMPENDING PLY FAILURE AFTER MULT IPL ICSTR06190 



SAT ION M. THE LOAD FACTOR:* I STR06200 

FORMAT <//lHO»2X, • CYCLE=* , I 3» 5X , » IFAI L=* t I3,5X t 'LOAD FACTOR= ' , G 10. 4STR06210 
l) STR 06220 

FORMAT ( 1H0»3X, ‘PLY* ,8X, ' STRESX ' , 8X, • STRESY* , 10X, ' JAUXY' ) STR06230 

FORMAT(4X,I3,6X,G10.4,4X,G10..4,4X,G10.4) STR 06240 

FORMAT( /////// lH0f2X»' *** THE STRENGTH RATIO FOR THE LAMINATE IS'STR06250 



StG10.4) 
END 



STR06260 

STR06270 

STR06280 



/ 
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Appendix C 



Instructions for Use of the Computer Program 



The following describes the sequence of input data for the 
computer program: 

1. Read (5,40) NPLY , I TYPE, ZNUXY, EX, EY, ALFAX, ALFAY, 

STFRTM 

40 Format (213, 7G10.4) 



2 . 



3. 



4. 



NPLY = Number of plies above the laminate midplane 

ITYPE is 0 if all plies are of equal thickness; 

ITYPE is 1 if plies have different thicknesses 

ZNUXY (v ) = lamina longitudinal Poissons ratio 
xy 

EX ( E x ) = lamina longitudinal stiffness 
EY (Ey) = lamina transverse stiffness 
ALFAX (a^) = lamina longitudinal coefficient of 
thermal expansion 

ALFAY (a ) = lamina transverse coefficient of thermal 
Y J 

expansion 

STFRTM = stress free temperature (i.e., the curing 
temperature of the laminate) 

Read (5,110) THICK (not inputted if ITYPE = 1) 

110 Format (G10.4) 

THICK = laminate thickness 

Read (5,120) (J,PLYTHK (I), I = 1, NPLY) (not inputted 
if ITYPE = 0) 

120 FORMAT (13, G10.4) 

J = ply number 

PLYTHK (I) = ply thickness 

Read (5,130) (J, ANGLE (I), TEM (I), I = 1, NPLY) 

130 Format (13, 2G10.4) 



J = ply number 

ANGLE (I) = <{k = ply orientation angle 
TEM (I) = T i = Ply temperature 
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5. Read (5,140) ZN1, ZN2 , ZN3 
140 Format (3610.4) 

ZN1 = = normal force along laminate 1-axis 

ZN2 = N£ = normal force along laminate 2-axis 
ZN3 = = shear force on laminate 

6. Read (6,140) ZM1 , ZM2 , ZM3 

ZM1 = M^ = bending moment on laminate 1-face 
ZM2 = M 2 = bending moment on laminate 2- face 
ZM3 = = twisting moment on laminate 

7. Read (5,170) X, XPRIME, Y, YPRIME, SS, FXYSTR 
170 Format (6G10.4) 

X = longitudinal tensile strength of lamina 

XPRIME (X*) = longitudinal compressive strength of lamina 

Y = transverse tensile strength of lamina 

YPRIME = (Y') = transverse compressive strength of lamina 

SS (S) = lamina shear strength 

FXYSTR (F* ) = O.S (in accordance with Eq. (59)) 
xy 

A sample output of a typical analysis is presented on the 
pages following this page. This is the problem of a cross-ply 
E-glass/epoxy laminate subjected to uniaxial loading, i.e. 

(N^, N 2 , N^) = (1,0,0). The cross-ply ratio, M, is 0.2. The 
curing temperature is 270°F, and the temperature of the lamin- 
ate is 70°F. 
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STRESSES AT IMPENDING PLY FAILURE AFTER MULTIPLICATION W. THE LOAD FACTOR t 

PLY STRESX STRESY TAUXY 

2 . 15000+06 -. 26580-09 .0 
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